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Abstract 

We develop an iterative, numerically exact approach for the treatment of nonequilibrium quan- 
tum transport and dissipation problems that avoids the real-time sign problem associated with 
standard Monte Carlo techniques. The method requires a well-defined decorrelation time of the 
non-local influence functional for proper convergence to the exact limit. Since finite decorrelation 
times may arise either from temperature or from a voltage drop at zero temperature, the approach 
is well suited for the description of the real-time dynamics of single-molecule devices and quantum 
dots driven to a steady-state via interaction with two or more electron leads. We numerically 
investigate two non-trivial models: the evolution of the nonequilibrium population of a two-level 
system coupled to two electronic reservoirs, and quantum transport in the nonequilibrium Ander- 
son model. For the latter case, two distinct formulations are described. Results are compared to 
those obtained by other techniques. 

PACS numbers: 03.65.Yz, 05.60. Gg, 72.10.Fk, 73.63.-b 
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I. INTRODUCTION 



There are several ways in which a quantal entity may exhibit nontrivial departures from 
equilibrium. First, a system may evolve toward equilibrium after application of a transient 
external pulse or from a nonequilibrium initial condition. Simple examples of such situations 
are now relatively well understood [l, 2]. Related to these types of departures from equilib- 
rium, but less well understood, are more challenging cases of "quantum quenches," whereby 
the sudden change of a control parameter induces dynamics that probe non-trivial aspects 
of strong correlation or quantum criticality [3]. Also underdeveloped is our understanding of 
quantum mechanical systems driven to nonequilibrium steady-states via coupling to two or 
more electronic reservoirs. Since this is the case of direct relevance for the study of transport 
through quantum dots and molecular electronic devices ^,5], the complete description of 
this type of nonequilibrium behavior is of practical as well as fundamental interest. 

There are essentially two main theoretical frameworks for the calculation of properties 
related to the approach to, and attainment of, nonequilibrium steady-states of the types 
mentioned above. The first is the standard real-time Schwinger-Keldysh technique (| . This 
approach has led to the exact formulation of steady-state properties (e.g. the current) in 
terms of Keldysh Green's functions [7]. A variety of direct perturbative and renormalization 
group calculations have naturally emerged from this starting point 8NlO||. In addition, real- 
time Monte Carlo methods have been formulated on the basis of the Schwinger-Keldysh 



approach 
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The Monte Carlo methods are exact in principle but may be severely 
limited by numerical sign problems, depending on the formulation, system and regime under 
investigation. 



The second framework involves the use of Lippmann-Schwinger scattering states 



ig to 



construct the properties of nonequilibrium quantum steady-state. This approach has led 
to several rigorous results for integrable models 17j . In the last few years this viewpoint, 
combined with the notion of Hershfield's steady-state density operator 18|, has inspired 



the formulation of new non-perturbative approaches as well as numerical methods [19M22] . 
Most recently, promising numerical renormalization group approaches have been put forward 
based directly on the construction of scattering states 23), and an extension to the density 
matrix renormalization group method, incorporating real-time evolution, has been presented 



24. 



25]- 
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Consideration of more standard classes of nonequilibrium relaxation in dissipative sys- 
tems such as the spin-boson model has led to a variety of path-integral techniques for the 
numerically exact propagation of the reduced density matrix of a small system cou pled to its 
environment [l| . These methods, which include real-time Monte Carlo techniques llMl3| as 



well as deterministic iterative approaches 26H28] . are connected to the Schwinger-Keldysh 



type framework discussed above. Here, as in the Schwinger-Keldysh technique, the approach 
to equilibrium along a particular time contour from a prescribed nonequilibrium initial con- 
dition is described. Of these approaches, iterative path-integral methods have had particular 



success 



26j . Such methods are based on the notion that a well defined bath correlation time 



(if one exists) renders the range of the influence functional (IF) finite, allowing for a con- 
trolled truncation of memory effects and thus a deterministic propagation of observables 
that is free of the real-time sign problem. 

While iterative path-integral approaches have been proven successful in describing 
nonequilibrium dynamics in simple spin-boson type models in the last 15 years, only recently 
they have been formulated and used in cases of relevance to transport through quantum dots 
and molecular electronic devices 29]. In such systems, given that a chemical potential dif- 
ference between electronic reservoirs leads to a well defined decorrelation time for dynamics 
even at zero temperature, a memory time, beyond which correlations can be dropped, exists. 
This finite-memory characteristic allows the development of iterative techniques, capable of 
describing relaxation in a wide, non-trivial region of parameter space. 

In this work we develop and apply a new iterative path-integral technique to two models 
of nonequilibrium transport and dissipation: the spin-fermion model and the single-impurity 
Anderson model. The techniques developed here hold the potential for the exact description 
of long time dynamics in systems driven to a nonequilibrium steady-state via coupling to two 
or more electronic reservoirs. The method we describe in this work is conceptually similar 
to the ISPI approach of Thorwart, Egger and coworkers {29]]. The distinction between these 
two approaches lies mainly in the propagation scheme and the manner in which the leads 
are traced out of the problem. In the iterative approach developed here, the reservoirs 
are represented as discrete levels and are eliminated numerically via the Blankenbecler- 
Scalapino-Sugar (BSS) identity 30J. While this approach has the disadvantage that an 
additional source of systematic error is introduced due to the discretization of the lead 
degrees of freedom, we find empirically that the error is easily controlled without undue 
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computational expense. The advantage of this approach is that the study of general models 
(for example multi-site Hubbard "dots" ) may be performed with essentially no reformulation 
of methodology. Taking advantage of this fact, we present a first set of exact results for 
the out-of-equilibrium two-lead spin-fermion model. A second difference between the ISPI 
approach and the approach outlined in this work is related to the propagation scheme. Here 
we combine our matrix formulation with a propagation scheme similar to that described 



in 



This allows for very efficient propagation that may be trivially parallelized with 



commercially available software 



31] . These distinctions in scaling and flexibility of approach 



render our formulation as a useful compliment to the previously developed ISPI method. 

This paper is organized as follows. Section II and Appendix A present some general 
aspects of the iterative propagation technique. Section III contains a case study of the 
relaxation of a tunneling system coupled to two electronic reservoirs. In Section IV we in- 
vestigate nonequilibrium transport through an Anderson dot. In Section V we conclude. 
We include an alternative formulation of our approach for the nonequilibrium Anderson dot 
in Appendix B. This formulation may also hold promise in related path-integral approaches 
such as the ISPI approach. Appendix C describes extensions to finite temperatures. Fi- 
nally, Appendix D discusses some aspects of the convergence analysis which is necessary for 
elimination of the systematic errors in the method. 



II. GENERAL FORMULATION OF THE ITERATIVE APPROACH 

We consider a generic many-body system, consisting of a finite interacting region coupled 
to two infinite non-interacting reservoirs. The Hamiltonian H can be partitioned into a 
zeroth order term H whose solution can be exactly obtained, typically containing few-body 
interactions, and a higher order interaction term H\. We introduce our iterative approach 
using the reduced density matrix, p$ = Tre{p}, obtained by tracing the total density matrix 
p over the reservoir degrees of freedom. The time evolution of ps{t) is exactly given by 

Ps (s",s';t) = Tr B (s"\e- lHt p(0y Ht \s'). (1) 

We decompose the evolution operator into a product of N exponentials, e lHt = (e iHSt ) ; 
St = t/N, and define the discrete time evolution operator Q = e tH6t . Different Trotter de- 
compositions can be employed for splitting this operator. For example, we find it convenient 
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to approximate Q ~ e tHiSt/2 e iH st e iHiSt/2 w j ien studying the spin-fermion model (Section III), 
while for the Anderson model (Section IV) we find that it is useful to employ a decomposi- 
tion of the form Q ~ e iH st/2 e iH 1 st e iH 6t/2 ^ rp^g overg j] ^j me evolution can be represented in 
a path integral formulation, 

ps(s", s',t) = J dsQ J dsf... J ds%_ x J ds J ds±... J ds N _ 1 
^{(sVl^Xs+^iS^^ (2) 

where are subsystem (or fictitious) degrees of freedom, representing the discrete path on 
the forward (+) and backward (— ) contours. As an initial condition we may assume that 
p(0) = pbPs(O) with the bath (B) uncoupled to the subsystem. In what follows we refer to 
the integrand in ([2]) as an "Influence Functional" (IF) 
assigning sj^ = s", s^ = s'. Note that our definition o: 
contained in the original work of Feynman and Vernon .. 
to make connection with the iterative schemes developed in the previous path-integral based 



321 ]. and denote it by I(sJ, sf...s^), 
re IF is more general than that 
. We chose this loose definition 



numerical work 



26|. 



The IF combines the information of subsystem and bath degrees of freedom with system- 
bath interactions, and its form is analytically known only in special cases. For example, for 
a harmonic bath bilinearly coupled to a subsystem the IF is an exponential of a quadratic 
form, multiplied by free subsystem propagation terms 

N k 

i4-- s %) = ex P [- $^$^(4 - s k)(vk, k ' 4 -vl, k ' s k') 

k k'=0 

x (4|e-™l4-i)---(4lPs(0)i«o)---(^-i|e™|^}. (3) 



rhar I 



The coefficients 77^' depend on the bath spectral function and the temperature [26]. For 
a general anharmonic environment the IF may contain multiple-site interactions, where 
the coefficients are not known in general [271 ] . However, even when the form of the IF is 
analytically known as in ([3]), it still combines long range interactions limiting brute force 
direct numerical simulations to very short times. 

For a system coupled to a single thermal reservoir this challenge has been tackled at finite 
temperatures where a natural bath decoherence time exists. As noted by Makri and Makarov 



261 ]. such cases are characterized by the useful feature that nonlocal correlations contained 



in the IF decay exponentially, enabling a (controlled) truncation of the IF that includes only 



a finite memory length. Based on this feature, an iterative scheme for evaluating the (finite 



dimensional) path integral has been developed 26j . While the original quasi-adiabatic path 
integral (QUAPI) algorithm was developed based on the analytical pairwise form of the IF 
specific to harmonic reservoirs ([3]), a subsequent more general approach proposed in Ref. 



is based only on the fact that memory effects at finite temperatures generically vanish 
exponentially in the long time limit. 

This idea can be further employed to simulate the dynamics of a generic nonequilibrium 
bias-driven system [29|. Since in standard nonequilibrium situations bath correlations die 
exponentially, the IF can be truncated beyond a memory time r c = N s 5t, corresponding 
to the time where beyond which bath correlations may be controllably ignored. Here, N s 
is an integer, 5t is the discretized time step, and r c is a correlation time dictated by the 
nonequilibrium situation. For a system under a dc potential bias Afj, at zero temperature, 
r c ~ 1/A/i, while at temperatures for which T > A/z temperature sets the scale of the 
memory range. We therefore write the total influence functional approximately as 

with 

t( \_ I( S k ' S k+V S k+N s ) / K \ 
l s \Sk, Sk+NJ — ± ± y [0) 

1 \ S k i S k+li ••■■> S k+N s -l) 

The errors in Eq. (j3J) are the usual Trotter error arising from the time discretization and the 
truncation to a finite memory time r c = N s 5t. Both of these errors can be controlled. Eq. 
(j4]) can be understood as a simple generalization of the pairwise expression ([3]) for which 

i! ar (4,4 + v,4 + N s ) = fo(4)M4,4 + i)---fNs(4,4 +Ns )- (6) 

The one-body and two-body functions / can be obtained by rearranging Eq. ([3]). From 
these expressions we recursively build the finite-range IF for a general model. We assume 
that the complete functional decays to zero with time constant r c = N s 5t, (N s < N), thus 
it can be approximated by the product 

V i 1 ' 2 i '"' N—l) rf ± ± ± ± V v V 
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By recursively applying this rule, the truncated IF is further decomposed until it correlates 
interactions within r c only, 



I( s o ) S l! S 2 5- --5 S Ar)~ 

ji ± ± ± \ I( s l ; g 2 ; -1 S N s +l) H s 2 ) g 3 ) g Af s +2) H S N-N s i S N-N S +V ■■■■> 4) /q\ 

\1' 2'"'' N s ) 1 \ b 2 i b Z i ■i b N a +l) 1 \ b N-N s , b N-N 3 +l, b N-l) 

resulting in Eqs. ffl) and (^j. The physical content of this approach, which is similar to 
that described in [28], is outlined in Appendix A. The approach becomes exact as r c — > oo. 
Outside of the initial propagation step, I s (s , s 1 , s N ) = I(sq , s x , s^- ), we can identify 
the functions J s [Eq. fll])] as the ratio between two IFs where the numerator is calculated with 
an additional time step, Eq. fl5]). Next, based on the decomposition (@| we can iteratively 
integrate Eq. ([2]) by defining a multiple-time reduced density matrix p~s( s k, Sk+i, ■■, $k+N a -i)- 
Its initial value is given by p~s(4 , _ x ) = /, and it is evolution is dictated by 

p 5 (s± s%J = J ds±p s (s±, s^.J/^sJ, s±J, (9) 

with 

J s (4,...,s^J = Tr B {(s+J^ t l4 s -i)---(4l^ t l4)(4 Ip(°)I s oX s o l^l s r)--( s v s -il^l s v s )}- 

(10) 

A general propagation step involves integration over two (±) coordinates, 

Ps{4+v-> s k+N s ) = J d 4ps (s^,-, 4+^-1)^(4, -,4+n s ), (11) 

where the time-local (tk = k5t) reduced density matrix is obtained by summing over all 
intermediate states, 

Ps(tk) = J ds^_ 1 ...ds^_ Ns+1 p s (s^_ Ns+1 , sf). (12) 

The evolution at shorter times k < N s can be calculated in a numerically exact way. Before 
turning to specific models we would like to make the following comments regarding the 
above derivation, (i) The specific partitioning of the Hamiltonian into H and Hi depends 
on the model investigated. As we show below, H Q may include only the subsystem degrees 
of freedom (spin-fermion model), or it may be constructed involving all two-body terms 
(Anderson model), (ii) Obviously, the decomposition (J7J) is not unique, however, different 



schemes should lead to equivalent time evolution, and thus the partitioning is a matter of 
numerical convenience, (iii) The truncated IF (J s ) is not necessarily a time invariant. As we 
show below, in the spin-fermion model I s does not depend on time, thus in this case it needs 
to be evaluated only once during the propagation scheme. In contrast for the Anderson 
model standard use of the Hubbard- Stratonovich transformation leads to an IF expression 
that has to be updated at each time step. In Appendix C we outline an approach that 
does not make use of the Hubbard- Stratonovich transformation and thus produces a form 
on the IF of the Anderson model that is time- independent, (iv) The short-range function I s 
can be analytically evaluated in some special cases {3, 0|. For general reservoirs it may be 
evaluated numerically, by using finite size reservoirs as described in the next section, (v) The 
approach outlined here is not restricted to specific statistics of the leads (boson or fermion) 
and is solely based on the fact that at finite temperature and/or finite bias bath correlations 
exponentially decay at long time. Therefore, it can be used to treat finite temperature 
anharmonic bosonic environments [28] as well as nonequilibrium Fermi systems. 



III. DISSIPATION IN THE NONEQUILIBRIUM SPIN-FERMION MODEL 



A. Model 



As a first example, we consider the dynamics of a two-state system coupled to two 
fermionic leads maintained at different chemical potential values, the "spin-fermion model" 
(SF). This model has been considered in a series of recent papers 33M36l|. and serves as a sim- 
ple, albeit non-trivial, example exhibiting the generic behavior associated with the approach 
to a nonequilibrium steady-state. In particular, at zero temperature the chemical potential 
difference A/z sets the essential energy scale for dephasing as is expected generically in more 
complex models such as the nonequilibrium Kondo model [37J]. It should be noted, however, 

the connection between the model studied here and the nonequilibrium Kondo model 371 

f] 

is more tenuous then that between the tunneling center model in equilibrium [1] and the 



standard (equilibrium) Kondo model 



38 



. We take as our Hamiltonian 



H SF = H Q + Hi ; 
Hq = H S ] Hi = H B + H S b- 



(13) 
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The bath Hamiltonian Hb is taken to be that of two independent leads (a=L, R) charac- 
terized by (spinless) free-fermion statistics with different chemical potentials, namely 

H B = ^e k c\ k c a>k . (14) 

a,k 

The operator c a k (c a ,k) creates (annihilates) an electron with momentum k in the a-th lead. 
The system Hamiltonian Hs consists of a two-level system (TLS) with a bare tunneling 
amplitude A and a level splitting B, 

B A 

H s = -° z + ~o x . (15) 
We take the general form for the system-bath coupling to be 

HsB = Va,k;a' ,k'C a ,k C a' "*• (16) 



a,a'k,k' 



orms of the coupling parame- 



34 



391 ]. where the momentum 



Different versions of the model may be expressed via different 
ters V. In this paper we focus on the model presented in Ref. 
dependence of the scattering potential is neglected. The system-bath scattering potentials 
are then given by V a>a i, where a, a' = L,R are the Fermi sea indices. 

In the standard application of iterative path-integral approaches, two features greatly 
simplify the propagation algorithm. First, the form of the Feynman- Vernon influence func- 



tional is known analytically. Second, the influence functional is pair- wise decomposable 26 ]. 
As discussed in the previous section, neither of these features is necessary for the numerical 
implementation of an efficient iterative routine. 

Recently, the analytical structure of the influence functional in the spin-fermion model 
considered here has been elucidated, with a modified pair-wise Coulomb gas behavior emerg- 



ing at long times 35(. However, our recent numerical results have illustrated that in some 



cases for strong coupling of the system to the leads, most of the relevant dynamical evolution 



occurs in time intervals before strict Coulomb gas behavior holds [34]. 

The exact dynamics follows Eq. (T5]). Assuming separable initial conditions p(t = 0) = 
ps(t = 0)ps(i = 0), we can identify the IF in the present model as 

1(4, sf,..., 4) = (s+\p s {Q)\so)K{ S %,s%_ 1 )...K{4,8t)K(8t,8t) x 

Tr B | e~ iHl (4) 5 */2 e -i#i . . e -i#i («o )<W 2 pB (Q)e iHl ( s o")< 5 */ 2 .... e ^i(^-i)< 5 * e ^i(^)< 5 */ 2 j(i7) 
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where Hi = Hb + H$b provides an adiabatic partitioning of the Hamiltonian, s k 
are forward (+) and backward (— ) spin states along the paths, and K(s k+1 ,s k ) = 
{ s k+i\ e ~ lHsSt \ s t)( s k\ elHsSt \ s k+i) * s the propagator matrix for the isolated subsystem. 

The reduced density matrix is time-propagated by employing the iterative scheme flUJ)- 
( lT2"j) . where the function I s [Eq. (J5])] is calculated by taking ratios of the corresponding 
truncated IF (fTTj) . Note that this function is time-translationally invariant, thus we need to 
calculate it only once. 



B. Results 



To numerically calculate the influence functional, we express the lead Hamiltonians in 
terms of a finite number of fermions. Then, as in the standard BSS Monte Carlo approach to 
lattice fermions [30], the resulting trace may be expressed as a simple determinant containing 
the 1-body matrices that represent exponentials of operators that are quadratic in fermionic 
creation and annihilation operators. It should be noted that this discretization of the bath 
leads to systematic error in the results, unlike the case for the related ISPI approach of 



Thorwart, Egger and coworkers 



29j . However, the discretized approach for tracing out the 



bath is more flexible in that cases where the analytic structure of the self-energy terms, 
such as structured "dot" with several correlated sites, may be easily treated. Furthermore, 
bosonic analogs of generalized Anderson models may be treated easily as well 40], using the 



boson version of the BSS formalism 
developed bosonic versions of DMFT 



41 L This fact may be of importance for the recently 



42 



431 ]. where for out-of-equilibrium situations or at 
finite temperatures the approach outlined here may potentially serve as a real-time impurity 
solver. Fortunately, since the time intervals over which the bath is "measured" are short, we 
have found that the infinite bath result is easy reached even with a relatively small number 
of effective bath fermions ~ 40. 

We use the following parameters: A = 1, B — 0, A/x ~ 0.5 — 2, and pV aja > = A(l — 5 a ,a'), 
considering only inter-bath system-bath couplings, where spin polarization is coupled to 
scattering events between the nonequilibrium reservoirs. For simplicity we assume zero 
temperature. The generalization to finite temperature is straightforward as outlined in 
Appendix B. Since the iterative approach outlined above requires a finite range of memory 
for the influence functional, we work with a bias large enough to ensure facile convergence 
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in the numerical examples outlined below. 

In Fig. [1] we show the dynamics of the spin polarization (<J z (t)) for several different 
values of the bias A/i, distributed symmetrically between the L and R leads. The role 
of the chemical potential difference as a temperature-like contributor to dephasing is clear 



34(. We analyze (inset) the memory error in our algorithm by increasing r c , keeping St 



fixed. As expected, we find that r c roughly corresponds to 1/A/i. Thus, for A/i ~ 1, taking 
5t = 0.25, the dynamics is converging for N s > 5. A complete discussion of the appropriate 
convergence analysis is presented in Appendix D for the Anderson model. 




An=0.6 

---An=1.4 

Au-2.0 



30 



FIG. 1: Polarization in the nonequilibrium spin-fermion model at different values of the bias voltage 
A^=0.6 (full); A//=1.4 (dashed); A/i=2 (dashed-dotted) , B = 0, A=l, A=0.2, 5t = 0.25, N s = 8. 
Inset: convergence with increasing correlation time at A/i = 0.6, N s = 3 (dark full); N s = 4 
(dashed-dotted); N s = 7 (dashed); N s = 8 (dotted); N s = 9 (light full). Data was generated using 
80 states per bath, which is sufficient to ensure convergence in the regime of parameters presented 
here. 

In Fig. [2] we compare our numerically exact results, with the results of a generalized 
"non-interacting blip" approximation as formulated by Mitra and Millis 



33 



3J. While at 

weak coupling the dynamics reasonably agree, for strong interactions A = 0.3 (7cpV a ^ a ' ~ 1) 
the perturbative method diverges 34J. We found that at weak to intermediate interaction 
strengths our results systematically converge with increasing memory time r c . For strong 
interactions it 5 ~ 1 the time-step in our simulations should be made further smaller 5t ~ 0.1 
in order to achieve convergence, demanding extensive computation effort as N s > 16 for 
A/i ~ 0.6. 
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FIG. 2: Polarization in the nonequilibrium spin-fermion model at different spin-bath couplings, 
A=0.1 (dashed); A = 0.2 (full); A = 0.3 (dotted). Here A/j, = 0.6 and St = 0.25, N s = 10. 
The dotted line was generated using a nonequilibrium version of the "non-interacting spin-blip 



approximation" 



33, 



3j. 



It would be most useful to undertake a systematic study of the dynamical phase diagram 
in (T, A/i, pVay) space in the regions where our iterative technique is convergent. Such a 
study would be quite useful for the understanding of the approach to nonequilibrium steady 
state, and will be the subject of a future investigation. 

IV. NONEQUILIBRIUM TRANSPORT THROUGH AN ANDERSON DOT 



A. Method 



The single impurity Anderson Model (SIAM) 44| is one of the most important and well- 
studied models in condensed matter physics. While it was originally introduced to describe 
the behavior of magnetic impurities in non-magnetic hosts [38(, it has more recently served as 
a general model for understanding transport in correlated nanoscale systems j^, 5]. In such 
cases, the impurity is hybridized with more than one reservoir, and if the chemical potentials 
of the reservoirs are not identical, nonequilibrium transport will occur. Here, we present a 
numerically exact scheme for calculating dynamical quantities such as the time-dependent 
occupation and current in such systems. The approach outlined in this section relies on the 
discrete Hubbard- Stratonovich transformation. An alternative and more general approach 
is outlined in Appendix C that does not employ this transformation. While the approach 
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of Appendix C offers several advantages, it is somewhat simpler to implement the scheme 
described here, and for that reason we follow it for the sake of illustrative calculation. 

The SIAM model includes a resonant level of energy e d , described by the creation operator 
d\ (a =t, 4 denotes the spin orientation) coupled to two fermionic leads (a = L, R) of 
different chemical potentials p a , 

h am = J2 6d dl d a + Ud\ d t dld x 

+ Y £ k C l,k^,a + Y V ^kci A Ja + h.C. (18) 
a,k,a a,k,a 

Here k a (c a ^,a) denotes the creation (annihilation) of an electron with momentum k and 
spin a in the a lead, U stands for the onsite repulsion energy, and V a ^ are the impurity-a 
lead coupling elements. The Hamiltonian ( I18p can be also rewritten as H AM = H Q + Hi, 
where H includes the exactly solvable non-interacting part, and Hi includes the many-body 
term, 

(7 

+ Y e k C a,k,a C cc,k,cj + Y V ^ C l,k,a d ^ + h.C. 

Hi = U[n d ^n diX - -(n d|t + n d:i )] . (19) 

Here n d ^ = d\d a is the impurity occupation number operator. The shifted single- 
particle energies are denoted by E d = e d + U/2. We also define T = ^2 a T a , where 
r a = K Ylk \Va,k\ 2 5{ € — e k) is the hybridization energy of the resonant level with the a 
metal. 

Our objective is to calculate the dynamics of a quadratic operator A, either given by 
system or bath degrees of freedom. This can be generally done by studying the Heisenberg 
equation of motion of the exponential operator e XA with A here a variable that is taken to 
vanish at the end of the calculation, 

(A(t)} = Tr(pi) = 

hm ^-Tr \p(0)e iHAMt e xA e- iHAMt ] . (20) 

Here p is the total density matrix. For simplicity, we assume that at the initial time (t = 0) 
the dot and the bath are decoupled, the impurity site is empty, and the bath is prepared in 
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a nonequilibrium (biased) zero temperature state. The time evolution of A can be obtained 
following a scheme analogous to that outlined in Section II for the reduced density matrix. 
For clarity, we re-derive an explicit expression for the generalized IF in the present case as 
well. 

First we use a standard factorization of the time evolution operator e lH Mf = (e lRAM <5f ) iV , 
and assume the Trotter decomposition e lHAu5t ~ ( e %H oSt/2 e %HiSt e iH st/2\ man y Doc iy 



term Hi is further eliminated 
Stratonovich transformation 



introducing auxiliary Ising variables s = ± via the Hubbard- 



45| 



iHiSt _ 1 y ^ e ~sK+(n dtt -n dA ). 
s 

_ g-*W-(ttd,t- n »t,4-) ) (21) 



where k± = k' =f ik", k' = sinh 1 [sm(6tU/2)] 1 / 2 , k" = sin 1 [sm(StU/2)} 1 ^ 2 . The uniqueness 
of this transformation requires USt < n. In what follows we use the following notation 

e H±(s) ^ e -sK±{n dA -n dti ) ^ ^22) 

Incorporating the Trotter decomposition and the HF transformation ([221) into Eq. (|20|) . we 
find that at zero temperature the time evolution of A is given by 

(Mt)) = 

Jim — (0| ( e iH ^t/2 e iH 1 St e iH 5t/2\ N e \A / e -iHo<5t/2 e -iH 1 «5t e ~i/fo^/2\ Ar |q\ 

A^O <9A ^ 1 ^ ' ^ / I / 

= lim A|J_ y ds ± rfs ± ..ds±(0\ ^Ho/2St e H + (s+) e iH /23t^ ^iHoSt/2 e H + (sl) jHoSt/2^ 

where |0) is the initial (zero temperature) state of the total system. For convenience, we 
evaluate Eq. ( 1231) by diagonalizing the Hamiltonian Hq [see Eq. ( 1191) ]. and rewriting H± in 
terms of the new basis 

H = ^2e v b%; Hv = VH a V~ l 

H± = X>^' &t A', (24) 

with /5„ as the transformation matrix elements. We further transform both the operator of 
interest and the ground state into the new representation A = VAV -1 , |0) = V -1 |0). The 
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IF is identified as the integrand in (|23|) . where we truncate interactions beyond the memory 
time r c = N s St, 

I( s k ' ••■■ s fc+7v s ) = 

(25) 

with G+(s^) = (eiUo6t/2 e H + ( s ±) e iH 5t/2^ and g _ Finally, we can build the function 
I s [Eq. using fl5]), and the operator of interest A may be propagated using a scheme 
analogous to that developed for the reduced density matrix, Eqs. (!9|)- (fT2l) . 

Before presenting numerical results we make the following comments. First, in the present 
scheme the IF needs to be updated at each time step since the truncated IF [Eq. (1231) ] 
explicitly depends on the present time tk = k5t. Second, the operator A can represent 
various quadratic operators. Thus quantities such as the impurity population or the current 
through the junction [ijj may be investigated on the same footing. 

B. Results 

The IF (1251) is the core of our calculation. It is evaluated numerically using the zero 
temperature relationship (0|e B |0) = det[e fc ] occ ., where 6 is a single particle operator, B = ^ b, 
and the determinant is carried over occupied states only. Extensions to finite temperature are 
standard, see Appendix B. Similarly to the spin-fermion model we represent the reservoirs 
by a finite set of fermions, with energies determined by the metals' dispersion relation. 
Calculations must be converged with respect to the number of discrete lead states. The A 
derivative in ( 1231) is handled numerically, by calculating the IF for several (small) values of 
A. 

In the following we typically use the following conventions and parameters: a symmetri- 
cally distributed voltage bias between two leads with A/i = 0.4 — 0.6, a reservoir bandwidth 
of D = 1, a resonant level energy = 0.3, and hybridization strength T a =0. 025-0.1. 
Note that the actual hybridization parameter utilized in the simulations is the coupling 
V a ,k = \f^a/^Pa) where p a is the density of states of the a lead. For these parameters we 
find that convergence is achieved using L < 240 states per spin per bath. We have also 
verified that for A/i = 0.4 the memory time r c ~ 3.2 leads to convergence with 5t = 0.8 and 
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N s = 4, provided ^ < 3 (see Appendix D). 




25 30 



FIG. 3: Resonant level dynamics at different values of the voltage bias, A[i = 0.6 (dashed); 
A/jl = 0.4 (full). U = 0.1, T a =0.025, E d = 0.3, r c = 3.2. The dotted lines show for reference the 
exact U = dynamics at Afi = 0.6,0.4, (top to bottom). The circles are the respective Monte 
Carlo points. Calculations are performed at T = 0, while Monte Carlo data utilizes T = 1/200 
which is effectively converged to the T = limit. 

We begin by investigating the dynamics for a relatively small interaction U = 0.1 (r = 
+ r# and U/Y = 2). In this regime we are able to systematically converge the results of 
our procedure with respect to the three sources of systematic error, namely those associated 
with time step and bath discretization as well as non-local memory truncation. Figure [3] 
presents the time evolution of the dot occupation for two different bias voltages, A/x = 0.6 
(dashed) and A/j = 0.4 (full), assuming the dot (E d = 0.3) is initially empty. The results 
are compared to exact real-time Monte Carlo (MC) simulations employing the hybridization 



expansion 



46j manifesting good agreement at this relatively small U : At short times the 



IF data reproduce the MC features, while close to steady-state the MC results become 
increasingly unstable. The more recently developed weak-coupling expansion 47] is capable 
of significantly extending the time regime for which converged results may be obtained via 
Monte Carlo for symmetric cases, however this restriction limits the cases for which long-time 
results may be obtained. The MC data presented in this paper was generated at finite-low 
temperature, 1/T — 200. We have verified (data not included) that for this temperature 
range the population dynamics essentially coincide with the strictly zero temperature case. 
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The extremely small deviations between MC data and our approach at U = 0.1 in Fig. [3] 
are the result of small differences in temperature and the fact that a sharp, finite band is 
assumed in our calculations. 




5 10 15 20 25 30 
t 



FIG. 4: Population of the resonant level in the Anderson model. The results for U = (full), 
U = 0.1 (dashed), U = 0.3 (dashed-dotted), U = 0.5 (dotted) are compared with the exact 
dynamics at U = (o) and Monte Carlo data (*, □, <). The physical parameters of the model 
are D = 1, Afi = 0.4, Ed = 0.3 and r a =0.025. The numerical parameters used are L = 240 lead 
states, r c = 3.2 with N s = 4 and 8t = 0.8. Note that convergence and thus agreement with Monte 
Carlo cannot be achieved for t > 10 if S > 3. 



Figure S] presents the time evolution of (ri<i )0 -) with increasing on-site interaction. While 
we have not been able to overcome convergence issues for all times and all values of ^, 
we find that dynamics are faithfully reproduced for all jr at short times, while accurate 
and converged results are correctly obtainable only for ^ < 3. The strict requirements 
for convergence are presented in Appendix D. While this regime is one where perturbation 
theory in U is accurate 47J, |48| , we believe that convergence restrictions are surmountable 
within the methodology presented in this work. Future study will be devoted to this issue. 
Fig. compares the early propagation obtained within the IF approach (□) to the MC data 
(o). Interestingly, while our approach does not capture the t 2 characteristic at < t < 3 due 
to the rough time discretization, the intermediate time dynamics is still correct. It should 
be possible to devise an adaptive time propagation scheme where the time step is increasing 
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with time, keeping r c fixed. Future work will be devoted to improving convergence for large 
U and t. It is interesting to note that even though the results at large time and on-site energy 
(U/T > 3) are not converged and thus do not controllably represent a reliable estimate of 
population dynamics, the results are still reasonably close to the MC data even for ^ = 6. 




t 

FIG. 5: Short time dynamics in the Anderson model (□) compared with Monte Carlo data (o) 
for U = 0.5,0.3,0.1 (top to bottom). D = 1, A/x = 0.4, E d = 0.3, r Q =0.025, L=240, r c = 3.2 with 
N s = 4 and 5t = 0.8. 



V. CONCLUSIONS 

We have presented here a general path-integral based iterative scheme for studying the 
dissipative dynamics of bias-driven nonequilibrium systems. Our method relies on the finite 
range of bath correlations in out-of-equilibrium cases, thus interactions within the influence 
functional may be truncated beyond a memory time dictated by the nonequilibrium con- 
ditions, and an iterative and deterministic scheme may be developed. This scheme is in 
principle exact for cases where convergence with respect to truncation of memory effects is 
achieved. 

The philosophy of our approach is similar to the previously developed ISPI approach of 



Thorwart, Egger and coworkers 29j |. The distinction between the method presented here 



and ISPI is confined to the propagation scheme and the technique via which the leads are 



eliminated. The discretized BSS-like approach 30[ to tracing out the reservoirs used here 
may be employed in situations where the structure of the memory term is difficult to obtain 
analytically. Furthermore, the matrices involved in the iterative scheme are fixed in size, 
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and this fact may present numerical advantages at very long times. While our approach 
introduces an additional source of systematic error related to discretizing the leads, we have 
found that this error is easily controlled with limited numerical cost. Thus, our approach 
presents a related but complimentary methodology to the ISPI technique. It should be noted 
that currently the approach presented here and the ISPI technique appear to have difficulty 
converging in similar regions of parameter space that are accessible in some cases by, for 



example, the weak-coupling Monte Carlo approach [471 ] . However approaches like ISPI and 
the methodology presented allow for an accurate description of long-time dynamical features 
when they do converge, something that is generically difficult with Monte Carlo schemes. 
In this regard our approach is also complimentary to, and not competitive with, expansion 



based Monte Carlo schemes 



46 



47]. 



We have applied our technique to two prototype models: (i) The spin-fermion model of 
a spin coupled via a dipole-type interaction to two leads under a potential bias, and (ii) the 
Anderson model, where a resonant level with an onsite repulsion is coupled to nonequilibrium 
leads. In the first case the dynamics of the tunneling system was investigated, recovering 
damped oscillations for weak-intermediate couplings with the bias playing a role analogous to 
that of the temperature in equilibrium systems. For the nonequilibrium Anderson model we 
focused our study on the resonant level population. Our method yields results in reasonable 
agreement with numerically exact Monte Carlo simulations for weak to intermediate onsite 
interactions U. For strong U deviations are observed. The results presented in Appendix 
D suggest that the deviations are related to memory and time step truncation errors which 
we have been unable to control at the present time. Future work will be devoted to this 
issue. The study of more complex models, e.g. the multilevel Anderson model with onsite 
electron-phonon interactions will be the subject of future studies. 
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Appendix A: Justification of the Truncation Scheme 



Here, we justify the breakup of the IF as prescribed by Eq. (JSJ), demonstrating that 
the terms neglected account for interactions beyond the memory range r c . Consider for 
simplicity the functional 



, Si , So , So , Sa , Si 1 ~~ 



I( s ()i S l , S 2 i S 3 



I( s l , s 2 i s 3 > S 4 ) -^( s 2 : ■'"']■ -V, 



,S 



2 > ^3 



I{ S 2 > S 3 > S 4 



(Al) 
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with iV s =3. Using a cumulant expansion for the total 



'0 ' °1 ' °2 > °3 ' Cl 4 ' °5 

truncated here by following 
influence functional (IF) 
/ = 1(2) x 1(3) x 1(A) x 1(5), where each term is an exponent of a sum of the n-body 
terms, For example, 1(2) ~ e~^ i ^ 9i ' 3 with pairwise interactions g^j, 1(3) ~ e~^ i ^- k9l ' j ' k J 



28| , we write the IF as a product of n-body interaction terms, 



incorporating "three body" interactions gij,k- Substituting this structure into Eq. (lAljl . we 
find that the following terms are not present on the right hand side: The two- and three- 
body terms g 0A , #o,l,4, #0,2,4, #0,3,4, four-body terms, #0,1,2,4, #0,1,3,4 and #0,2,3,4, and a five-body 
element (70,1,2,3,4- These nonlocal interactions, connecting spins beyond the memory range 
specified, N s = 3, are assumed to be small, and are therefore discarded in our truncation 
scheme. Larger memory blocks, connecting more distant time slices, may systematically be 
included until convergence with truncation of memory terms is reached. 

To make this discussion concrete, consider a situation where non-equilibrium Coulomb 
gas behavior holds, as discussed in 



34 



35fl . In such cases, the total influence functional 



will be of the form / ~ exp 



J2i> j ~ where C (t) oc A[i\t\ up to logarithmic 



corrections. Consider now Eq. ([8]). Clearly the leading term contains all interactions 
between "charges" separated by a distance in time that does not exceed \to — tjvj, namely 



I(s%,sf,...,s% ) ~ 



exp 



Terms of the form 



i(st^-,4 include 

l(s 1 ,s 2 ,...,s N ) 



only interactions between "charges" interacting over the time intervals \t n — t ^+1 1 where 
< n < N B+ i, without double counting terms already contained in J(s , s 2 , s N ). This 
procedure is then iteratively continued until the complete influence functional is constructed. 
The error accrued originates from the neglect of terms in the exponent of the order A fir 
where r = \t a — t&| and b — a > N s + 1. Thus, the procedure is rendered controlled and 
is expected to converge to the exact result as long as N s is made large enough. It should 
be noted that the approach outlined here is more general than this and is expected to hold 
at short times or very large couplings where Coulomb gas behavior may break down, as 



discussed in 



35] 
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Appendix B: Extensions of the IF Technique to Finite Temperatures 



We present here the natural extension of our approach to finite temperature. The core 
of our numerical calculation is the influence functional (IF), incorporating the Fermi sea 
degrees of freedom, e.g. Eq. (II 7p for the spin-fermion model or Eq. ( 12 5 p for the Anderson 
model. Assuming for simplicity a single Fermi sea, consider the following IF-like object 

C f = Tr B [e M -e M *p B ] , (Bl) 

where M\ and M2 are quadratic operators and ps = e~^ HB /TrB[c~^ HB ], Hb is the bath 
Tamiltonian, (fl4|) . This correlation function can be expressed by single-particle operators 

C / = det[J-/(c) + e mi e raa /(c)]. (B2) 

Here /(e) = [1 + e - ^ e-At '] -1 is the Fermi-Dirac distribution function, f3 is the inverse tem- 
perature, I is the unit operator, and mi and m>2 are single-particle operators corresponding 
to Mi and M2 respectively. This expression can be trivially extended to include more ex- 
ponential terms, e Ml e 2 ... • e Mjv , as necessary for the evaluation of the IF expression. For 
multiple-independent reservoirs, p B = Pl® Pr, the above relation can be generalized, 

C f = Tr L Tr R [e Ml e M *p L ®p R ] 

= det {[(I L - f L (e)) ® I R ] [(/ - f R (e)) ® I L ] + e^e™ [f L (e) ® I R ] [f R (e) g) I L }} . 

(B3) 

Here is the identity matrix for the a space; a = L,R, and f a (e) = [1 + e"^"^)]- 1 . The 
above expressions reduce to the ones used in the text for T = 0. 



Appendix C: An Alternative Formulation: Nonequilibrium Transport Through an 
Anderson Dot 

We present here an alternative formulation for calculating the dot properties in the single 
impurity Anderson model (SIAM) without invoking the Hubbard-Stratonovich transforma- 
tion. This formulation is based on a different Trotter decomposition than that used in 
Section IV. While the resulting expressions are more complex for the decomposition de- 
scribed here, it has the advantage that the resulting IF need not be updated each time step. 
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Furthermore, since fewer terms of the Hamiltonian are split in the Trotter decomposition, 
it is possible that larger time steps may be taken with the decomposition presented here. 
Further work investigating this approach, which is not confined to the Anderson model, will 
be presented in a future work. We refer to the approach developed in Section IV as SIAM 
I, and to the method of this appendix as SIAM II. 

We begin by partitioning the Hamiltonian ( !T8|) as follows: Ho includes the subsystem 
(dot) terms, and Hi includes the two non-interacting leads {Hb) and system-bath couplings 
(H S b) 

H AM = H + Hi, Hi = H B + H S b 
H = ^ e d n djCT + Un^n^, 

<7 

H B = e kcl Aa c a:k , ai H SB = V ^cl Aa da + h.c. (CI) 

Here n d (T = d\d a is the impurity number operator and k a is a creation operator of an 
electron at the a lead with a spin a and momentum k. Note that H can be explicitly 
described by a 4-state system, |1) = |0,0), |2) = | f, 0), |3) = | 4,0), |4) = | t,l), corre- 
sponding to an empty dot, a single occupied dot of a =J [,i, an d a double occupancy state. 
When U is very large {U — > oo), we effectively have a 3-state system, since double occupancy 
becomes negligible. The energies of these four subsystem states are E\ = 0, E 2t 3 = e d , and 
E 4 = e d + U. 

Consider the reduced density matrix p s = Tr B {p} obtained by tracing the total density 
matrix p over the reservoir degrees of freedom. The time evolution of ps{t) is exactly given 
by 

p s (a,a',t) = Tr B (a|e- iHAM V(0)e iHAM *|a') ; (C2) 

where |a) and \a') are subsystem states, as described above. Using the standard Trotter 
breakup, e iHt = (e iHSt ) N , 5t = t/N, and e iH ' AMst w ^1/2^81^51/2^ we can rewrite Eq> 
12]) in a path integral formulation, 

ps(a,a\t) = J dsQ J dsf... J <is^_ 1 J ds$ J dsi...J c?s^ v _ 1 
Tr B {(a\e- iHoSt / 2 e- lH i St e- lHoSt / 2 \s+_J 

(sM0)\s }...(s N ^\e tHoSt / 2 e^ St e i ^ (C3) 
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where Sk are subsystem states. As an initial condition we may assume that p(0) = pePs(O) 
with the bath (B) uncoupled to the subsystem. We focus next on the following matrix 
elements in Eq. (1C3P 



G^St) = (a\e- iHoSt/2 e- iHlSt e- iHoSt/2 \b) = e-^ +E » )5t/2 {a\e~ iH ^\b) . (C4) 

To compute {a\e~ lHlSt \b) note that it is advantageous to use again the Trotter splitting 

(a|e-^ l5 *|6) ^e- iHBSt/2 (a\e- iHsB5t \b)e- iHBSt/2 . (C5) 

We thus focus next on the matrix element 

O a , b = (a\e- iHsB5t \b), (C6) 

a quadratic operator in the space of the non-interacting electrons. It is useful to define the 
"composite" fermion c 0i<7 = Y^ a ,k V a,kC a ,k,^ leading to H SB = (4,<A + <4 c o,a)- In this 
representation a direct expansion of the exponential gives 

e XHsB = I + (cosh A - l)a 2 + sinh Aai (C7) 

with A = —iSt, ai = Hsb and 0,2 = J2a {dad^c^^c^ad^d^ + h.c. j . The operator ( 1C6|) 
is therefore of the form, O a ,b = a + /3co jCT + + •yc^co^ + 7 / co j0 -Cq (T , with constant 

coefficients a,fl,j. Substituting the pieces (1C5j) - (1C7|) into Eq. flC4j) yields 

Ga,b(^) « e^a+W/V^*/ 2 ^-^"/ 2 , (C8) 

incorporating linear combinations of bath operators Cq 7<7 up to a quadratic order. Finally, 
we put all pieces together into Eq. (1Q3|) and obtain the reduced dynamics 



p s (a,a',t) = J ds£...j ds%_ x J ds ...J rf%_i(4lPs(0)|s ) 

N-l AT-1 

exp [~i8tJ2 E s j + & E % - ^ + ^4)^/2 + «(# a + ^-)^/2] 

3=1 3 3=1 

xTrJ (e- lHB5t / 2 O as+ e- lHB5t O s+ s+ e~ lHB5t ...O s+ s+e - iH * st / 2 ) p B (0) 

L-iH B 5t/2 -iH B 5t -iH B 5t q -W B 5t/2\ j (Cg) 

Identifying the integrand as the IF, we can use the approach of Section II, define the trun- 
cated IF I s , and iteratively propagate the reduced density matrix to long times. 



23 



The approach developed here (SIAM II) has three main advantages over the method 
described in the main text (SIAM I), see Section IV. First, since the present method does 
not rely on the Hubbard- Stratonovich transformation it can be applied to general many 
body interaction Hamiltonians, while SIAM I is restricted to the Anderson model. Second, 
since the Trotter error in SIAM II is due to system-bath factorization, rather than one- 
body- many-body splitting as in SIAM I, the method described here should be beneficial 
in calculating dynamics of weakly coupled system-bath models with arbitrarily large many 
body (local) interactions. Finally, this method also suggests a computational advantage over 
SIAM I, since the IF here [integrand of Eq. ( \C9\i ] is time independent, unlike the IF of Eq. 
(F2S]) which needs to be recalculated at each time step. 



imit is assumed and 

Q 



4(|. (ii) Trot- 



29. 



Appendix D: Convergence Analysis for the Anderson Model 

There are three separate sources of systematic error within our approach, (i) Bath dis- 
cretization error. The electronic reservoirs are explicitly included in our simulations, and 
we use bands extending from —D to D with a finite number of states per bath per spin 
(L). This is in contrast to standard approaches where a wide-band 
analytical expressions for the reservoirs Green's functions are adopted 
ter error. The time discretization error originates from the approximate factorization of the 
total Hamiltonian into the non-commuting H (two-body) and Hi (many-body) terms, see 
text after Eq. (120]) . While for U — > and for small time-steps St — > the decomposition is 
exactly satisfied, for large U one should go to a sufficiently small time-step in order to avoid 
significant error buildup, (iii) Memory error. Our approach assumes that bath correlations 
exponentially decay resulting from the nonequilibrium condition A/i / 0. Based on this 
crucial element, the influence functional may be truncated to include only a finite number 
of fictitious spins N s , where r c = N s 5t ~ 1/A/i. The total IF is retrieved by taking the limit 
N s A, (N = t/5t). 

These three errors can be systematically eliminated by increasing the number of bath 
states, choosing a small enough time-step, and adopting a sufficiently long memory time. 
Note however that the last two strategies are linked: Increasing r c essentially means in- 



creasing the time-step, since the memory length is restricted to sma 
for practical-computational reasons. Thus, as in standard QUAPI [26 



1 values N s = 4 — 6 
I, one should find an 
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optimal balance between the time-step error and the memory size that correctly represents 
the dynamics. Ref. jso| suggests a systematic approach for reaching convergence using the 
QUAPI method, eliminating the Trotter discretization error and the memory truncation 
inaccuracy by extrapolating the data to vanishing time-step and to infinite memory time. 

A similar idea can be adopted here. First the bath finite-size error can be eliminated 
by systematically increasing the number of fermions at each lead. As an example, Fig. [6] 
presents the dot population for U = 0.1 and U = 0.5 taking L= 20, 40, 80, 120 and 240 (top 
to bottom). The inset shows that convergence can be reached, and that the occupancy is 
systematically decreasing with L. Next, the Trotter error can be eliminated by extrapolating 
the data to the St — > limit. Fig. [7] presents as an example the occupancy for A/x = 0.4 
using r c = 3.2, and St = 1.6, 1.05, 0.8, 0.64. The inset manifests convergence as a function of 
(St) 2 . Note that in the asymptotic limit the data points are slightly enhanced, practically 
canceling the effect of the bath discretization. Finally, the memory effect is analyzed in 
Fig. El For the parameters employed here (Ed = 0.3, U = 0.1, T a = 0.025, A/i = 0.4) 
convergence is arrived at r c ~ 4 (inset), in agreement with the rough estimate r c ~ l/Afi. 
We have not been able to obtain full convergence for U/T > 3. 

Using this analysis, we have recalculated Fig. H] extrapolating our data to (i) L — > oo, 
(ii) St — > and (iii) r c — > oo. Since the extrapolations (i) and (ii), bring about counter con- 
tributions, see Figs. |6] and [TJ the overall effect of the bath-time step-memory extrapolations 
on the occupation is rather small, and Fig. H] remains essentially intact. 
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